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NUMERICAL  SOLUTION  OF  NONLINEAR  PARABOLIC  EQUATIONS 


Samuel  Schechter 


1 . Introduction 

We  will  consider  the  parabolic  differential  equation 


3g(8)  = 32f(9) 

at  ax2  (1.1) 


for  the  unknown  function  0(x,t)  over  the  set  0 < x < 1,  Oct^T 
and  given  boundary  and  initial  conditions 

9 (x,0)  = „c  (x)  (1.2) 

0 ( i 1 1)  = v.i  (t),  i = 0,1. 

We  assume  that 

g'(9)  > 0 (1.3) 

f '(9)  >0 

for  all  0 in  a set  S. 

Typical  examples  are  found  where 

m 

g(8)  = 8 , f (9)  = 0 , m > 0.  (1.4) 

For  m = 2 we  get  the  equation  for  flow  of  gas  in  a pipe  [ 9 ],  while 
Richtmyer  [ 4 J examines  the  case  m = 5 for  a running  wave.  Examples 


* 
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arise  in  flows  in  porous  media  [ 10].  Other  examples  may  be  found  in 
Ames  [ 12] . 

Wc  approximate  (1.1)  by  using  a weighted  form  of  Euler's  summation 
formula  (or  second  order  Pad6  approximation)  [13], 

g^  - g"  = tk(0!g”  + /3gj)  + (k2/12)(yg°  - 6g*)  (1.5) 


where 


J , J. 
g.  = g(u. ) , 


u.  = 0(x  ,t  + jk),  g = g 
l i t 


i = 1,2, ... , n;  j = 0,1 


nh  = l,k  = At,a  + j3>o  , 

and  <y,  ft,  y,  6 are  nonnegative  constants.  If  these  constants  are  specialized 

5 

to  all  equal  1,  then  we  obtain  the  usual  formula,  with  remainder  of  order  k . 

3 

Other  choices  give  a lower  order  accuracy,  usually  0(k  ). 

We  will  find  it  convenient  to  drop  the  superscript  1 in  later  formulas. 

To  obtain  an  approximation  for  (1.1),  we  replace  g by  f - in  (1.5),  where 

2 XX 

f - = (f  -2f  + f )/h  . This  will  yield  an  implicit  scheme,  with  a 

xx  1-1  i 1+1 

nonlinear  system  in  n unknowns  to  be  solved  at  each  time  step. 

If  we  take  the  special  case  (1.4)  and  set  Or  = y = 6 = 0,  £ = 2,  we  get 
the  backward  Euler  system 


u 


1 


0 
u 

1 


k(u 


i-1 


m , 2 
Ui+l/h 


(1.6) 


2 


SB** 


where  m is  here  an  exponent,  not  a superscript . 

To  simplify  our  notation  we  define,  for  the  vectors  u and  v in 

Rn , the  Schur  product  uv  = u o v = vu  as  the  vector  with  components 

u v . If  f(x)  is  a function  of  one  variable, we  define  the  vector  f(u)  by 
i i 

f(u)  = (f (u  ) f (u  ))T 

1 n 

also  called  a diagonal  mapping  [ 8 ]. 

In  particular, u"1  will  denote  the  vector  whose  components  are  raised 
to  the  power  m. 

For  each  vector  u we  may  correspond  a diagonal  matrix,  which  we 

denote  by  its  upper  case,  U = diag  (u  ,...,u  ) = diag  (u)  and  conversely. 

1 n 

It  then  follows  that  uv  = Uv, where  the  right  side  is  ordinary  matrix 
multiplication. 

If  f'(x)  is  the  derivative  of  f(x)  then  F,(u)  = diag  (f'(u))  is  the 
Jacobian  of  f(u).  If  u depends  on  t,  then  f(u)t  = f/(u)u  = F*(u)u. 

With  this  notation  in  hand, we  may  write  (1.6)  as 
cm2 

u - u = ~2\Ku  , 2X=k/h  , (1.6) 

where  K is  the  usual  tridiagonal  matrix  of  order  n 
K = [-1,2,  -1]. 

A1 ternatively,  we  may  write  (1.6)  in  quasilinear  form  as 
o 

u - u = -2XA(u)u  , 

. m-1 

where  A(u)  = KU 

3 


Returning  to  the  more  general  case  (1.1)  we  get 


g(u)  = -h  2 K f(u) 
g(u)  = -h  2 K F'  (u)u 

Since  g(u)  = G'(u)u, 
we  get,  after  solving  for  u, 

g(u)  = h L (u ) f (u ) , 


(1.7) 


where  we  set 


.-1.. 

L.  (U  / = H f \,UJ  U iuy  K 


Inserting  these  into  (1.5)  and  collecting  terms  ,we  get 

g(u)  + fyJKf (u)  + |X2 5 L(u)f(u)  = (1.8) 

o o,2  o o 

Rq  = g - Kf  + y h f , 

o 

where  the  zero  superscript  indicates  u is  replaced  by  u . We  seek  a 

solution  to  (1.8)  and  a method  to  construct  it. 

Since  the  latter  calculations  used  only  the  fact  that  K was  independent 
of  u and  t,  the  same  result  is  obtainable  for  a general  linear  operator  T 
so  that  g^  = T(f ) . 

We  now  assume  that  L(u)  is  replaced  by  Lq  = L(u)  where  u is  some 
approximation  to  u,*  a solution  of  (1.8).  If  we  denote  the  left  side 
of  (1.8)  by  p(u),  then  (1.8)  becomes 

r(u):  = p(u)  - RQ  = 0.  (1.8) 
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Lemma : Let  S = [0 : f ' (0  ) > o,  g' (9 ) > o] . If  S is  convex  (that  is  an 

* 

interval ), then  (1.8)  has  at  most  one  solution  u with  components  in  S. 

Proof : Assume  there  are  two  solutions, u and  v.  Then 

(u  - v,  p ( u ) - p ( v ) ) = o 

From  the  hypotheses  we  get 

(u  - v,  g (u ) - g ( v ) ) = (u  - v,  G'(z)(u-v))  > 0 

for  7.  in  the  interval  [u,v3,and  since  K is  positive  definite, 

(u-v,  K ( f (u ) - f ( v) ) = (u-v,  KF'(w)(u-v))  >0  . 

The  ipplies  to  L f.  This  gives  a contradiction  and  proves  the 

o 

L- 

we  assume  that  f(0)  and  g(0)  have  positive  derivatives  for  all  0 , 
then  we  obtain  the  existence  of  a solution.  From  the  previous  lemma  it 
will  be  unique.  Thus,  we  assume  that 

f'(0)  > 0,  g'(0)  > 0 for  all  real  0 . (1.9) 

We  may  then  assume  that  f(0)  s 0,  since  by  inverting  f(u)  = v 

we  get  a new  equation  in  v with  g(u)  = g(f  1(v).  With  f(u)  = u,  (1.8) 

becomes  a semilinear  problem,  with  g(u)  a diagonal  isotone  mapping. 

If  6 =0,  then  p(u)  becomes  an  M-function  [ 8 ] so  that  existence  is 

obtained  in  this  case  by  a theorem  of  Rheinboldt  [5  ].  In  the  case  where 

6 > o,  we  obtain  a gradient  map.  Since  the  Jacobian  of  p satisfies 

p'(u)  s 1|3K2  NSNnin00'1  ' 

5 


the  existence  and  uniqueness  follow  from  results  in  Schechter  [ 7 ]. 


Both  of  the  above  theorem  also  provide  methods  of  construction  of 
the  solution.  In  particular  we  may  use  a variety  of  relaxation 
methods  [8  J.  We  have  thus  proved: 

Theorem . If  f and  g satisfy  (1.9),  then  (1.8)  has  a unique  solution 
that  may  be  obtained  by  a relaxation  method. 

It  should  be  noted  that  the  solution  method  does  not  require  a 
priori  estimates  for  this  special  class  of  problems. 
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2.  Kxi stence  and  Estimates 

Although  existence  and  uniqueness  are  obtained  as  stated  in  the  previous 
section  the  actual  methods  available  to  compute  the  solution  are  usually 
dictated  by  the  nature  of  the  problem.  If  6=0, we  may  resort  to  nonlinear 
SOR  with  overrelaxation  parameter  lc  in  (0,1]  or  use  a Jacobi  iteration  [8  ]. 
Since  these  methods  tend,  in  general,  to  be  slowly  convergent  and  since 
we  wish  to  cover  the  more  general  case  6 > 0,  we  consider  the  use  of 
relaxation  or  SOR-Newton  methods  [6  ],  [7  ],  [14]. 

For  this  purpose  we  find  it  necessary  to  require  that  f and  g have 
continuous  derivations  and  at  times  continuous  second  derivatives.  This 
will  allow  us  to  choose  relaxation  parameters,  at  times,  outside  the  unit 
interval . 

Before  proceeding,  we  note  some  estimates  of  the  solution  for  the 
case  6=0.  We  may  assume  as  above  that  f is  the  identity  map  and  g(0)=0. 
Otherwise,  we  may  subtract  a constant  from  g. 

If  u*  is  the  solution,  it  then  follows  from  Taylor's  expansion  that 

(D  + A)u*  = Rq  (2.1) 

where  we  set  A = X£K  and  for  some  z in  [0,u*],  D = G'(z).  Since  D + A is 
a monotone  matrix,  its  inverse  is  nonnega ti ve  and  we  get 

u*  = (D  + A)-1  Rq  , (2.2) 

and  if  RQ  >0,  then  u*  >0. 

In  practice  we  may  achieve  RQ  > 0 by  choosing  either  a small  a or 
small  time  step,  if  g°  > 0.  in  this  case  we  also  obtain  an  upper  bound 
on  u*,  since  D > 0 , 

0 < u*  * A-1  Ro  (2.3) 
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We  note  in  passing  that  (2.3)  can  be  used  as  a basis  for  the  usual 
fixed  point  iteration  process.  Since  this  is  usually  slow  we  do  not 
pursue  it. 

A similar  upper  bound  is  available  if  6 S 0 since  L0  is  positive 
definite  and  monotone.  if  we  set  B = sA6l,,  then 

(D  + A + B)u*  = R (2.4) 

o 

and  if  u*  > O.we  get 

u*  < (A  + B)-1  R . 

o 

To  achieve  the  positivity  of  the  solution  we  may  choose  6 or  k small 
enough  so  that  RQ  - Bu*  > o.  We  may  of  course,  without  this  restriction, 
obtain  a bound  on  (u^.u*)  by  using  the  energy  function  associated  with 
the  symmetric  problem.  This  bound  yields 

|ro|  2 i ^min  (A)-  | u*| 

For  some  special  choices  of  and  6 where  6 is  not  small  this 
positivity  may  again  be  obtained.  If  ^3=2  and  6 = 3, then,  absorbing  X 
into  K,  the  operator  on  u*  may  be  written  as  D + 2K  + KLC}K  = (I  + KIto)(D  + K)  + 
K(I  - L0D) . Note  that  the  first  tern  on  the  ri^ht  is  monotone  since  the 
factors  are  M-matrices.  If  z and  u are  close, then  LQD  is  near  the  identity. 

If  the  time  step  is  small, this  will  prevail  with  smooth  functions.  If  the 
second  term  is  combined  with  R0,  an  argument  similar  to  the  previous  one 
may  be  applied  to  achieve  the  result. 

For  the  example  (1.4)  cited  above  we  see  that  the  derivative  of  f 
does  not  appear  to  be  positive  for  all  0.  In  many  applications.it  is 
usual  to  have  the  solution  nonnegative  so  that  if  f is  restricted  to  this 
set, we  do  have  a positive  derivative.  If  a solution  appears  in  that  set, 
we  showed  previously  that  it  will  be  unique. 
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In  order  to  achieve  existence  we  may  extend  the  definitions  of  f or  g 


so  that  the  positivity  conditions  are  valid.  If  a positive  solution  was 
available  before  the  extension  it  should  not  be  lost. 

Let  us  consider  the  example  m = 2 since  it  is  the  simplest  and  most 
common  case.  Other  cases  may  be  treated  in  a similar  fashion.  By  inverting 

X 

f we  get  g(9)  = 95.  We  now  extend  g by  defining 

g(S)  = | e|  sgn  (0  ) (2.5) 

This  function  has  positive  derivatives  except  at  9 = 0 where  it  is 
infinite.  The  function  is  -lonotone  so  that  uniqueness  is  available  as 
indicated  above.  To  get  existence  we  make  use  of  the  potential  function 
E(u)  of  r(u)  . This  potential  is  strictly  convex  and  E -♦  <*>  as  u -»  ®. 

Thus  its  level  sets  are  bounded,  yielding  a global  minimum  and  the 
existence  of  a solution. 

If  we  perturb  g by  a small  amount  we  can  get  a more  constructive 
method  of  solution  and  allow  for  a large  variety  of  solution  methods.  For 
example,  we  may  change  g to  be 

g(e,9)  = ( ( 1 9 1 + e)2  - e5)sgn  9 (2.6) 

How  does  this  change  the  solution?  We  assume,  in  general,  that  for 
all  9 > 0 , 

0 < g(  ex  ,9 ) - e(e2,9)  £ C(e2  - ep  for  < e2  (2.7) 

for  some  positive  constant  c depending  on  9.  Consider  first  the  case  of 
6=0  and  let  u and  v be  the  positive  solutions  to  the  equations  (1.8)  for 
the  parameters  and  e2,  respectively.  Thus, 

p(€j,u)s  g(eitu)  + Au  = g(e2,v)  + Av  s p(e2,v)  = Rq 

and 

P(ej , v)  > p(e2, v)  = p(eltu) 

9 


so  that  v > u.  This  follows  from  the  fact  that  p is  inverse  isotone. 

Since  the  solutions  are  bounded  from  below  by  zero  we  get  a decreasing 
sequence  of  solutions  as  e goes  to  zero.  This  again  verifies  the  existence. 

From  the  hypotheses  there  is  a positive  diagonal  matrix  C such  that 
C(e2  - €j)  ^ g(ej,v)  - g(e-,.v)  = g(ej.v)  - g(€j,u)  + A(v-u)  2-  A(v-u) 

A"lc<«2"el)  * V'U  > °‘ 

Thus,  if  we  let  ej  go  to  zero, 

A-1Ce2  2 v-u*  > 0 , 

which  estimates  the  perturbed  solution,  with  possibly  a new  positive  C. 

If  we  choose  e = o(km)  for  a suitable  m then  we  stay  within  the  error 
of  the  equation. 

Tlie  example  (2.6)  given  above  is  differentiable  for  all  6 and  satisfies 
(2.7). 

Consider  the  case  of  6 2 0 and  let  E(e,u)  be  the  potential  function 
of  r corresponding  to  g(e,u).  Thus,  r is  the  gradient  of  E,  We  assume  that  for 

alt  9 | B<V0)  - g(er6)l  * c|  e2  " eX  | • 

for  some  positive  c.  It  follows  that  the  integral  g of  g as  a function  of 
9 also  satisfies  such  an  estimate. 

Let  u^  and  u^  be  solutions  corresponding  to  ^ and  e ^ and  set 
E ( e ^ > u j ) = E . Then,  since  the  first  solution  minimizes  El^.u), 

E12  " Ell  = g(u2  ~ ui)T  e"<z><u2  - V 
2 t \nin(A)  I U2  - U1  I 2 ■ 
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We  will  show  that  the  left  side  is  of  order  ~ &1  = e'  First, 
Eij  " Eij  depends  only  on  g,  this  difference  is  0(e).  Since  E^. 

E12  " EH  = E12  “ E22  + E22  “ E21  + E21  “ Ell 

S E12  - E22  + E21  - En  ^ C le|  • 

This  yields  the  estimate 


If  e<2  -*  0>  it  follows  readily  from  the  uniqueness  that  u-j  -*  u . 
This  gives  the  estimate  desired, 

|u*  - uj  c 


since 


11 


3 . Quadratic  System 


We  consider  the  specific  example  (1.4)  with  m = 2.  We  then 
get  a quadratic  system  to  solve  for  given  Rq  > 0 

0 (v)  = v + Av2  - Ro  (3.1) 

with  A positive  definite  and  monotone.  The  Jacobian  of  p is 

p ' = I + 2AV  = M 

and  if  we  have  V 0,  then  p is  nonsingular.  It  will  be  soteven  if  some 
but  not  all  entries  in  V vanish.  Assume  that  v°  is  positive  and  set 

/ fX 

r (u  ) = M.  We  may  try  to  use  a modified  Newton  iteration  or  a block 
relaxation  method  to  find  the  positive  solution. 

The  iteration  appears  as 

M ( v1  - v°)  = -U>p(v°)  (3.2) 

for  a given  choice  of  parameter  We  assume  that  tc  ischosen  sufficiently 

small  so  that  v*  > 0.  If  o < 0,  then  since  M is  monotone  we  do  not  need  to 
restrict  t*.  in  this  case. 

Since  Mis  not  symmetric  0 is  not  a gradient  mapping  and  it  is  not 
apparent  that  this  iteration  will  converge.  If  we  transform  the  variable 

O 

to  u = v , then  a variety  of  choices  are  available  for  b.  to  yield  convergence. 

Let  us  examine  the  Newton  step  for  the  transformed  problem  for  u°  ^ 0 

r(u)  = + Au  - RQ  = P(v)  . (3.3) 

/ 

Let  J = r so  that 

J ( u° ) ( u 1 - u°)  = - W’r°  = - cppo 

J = tlT2  + A = £(I  + 2AV  )V“l  = SMV-1 
2 

I f we  replace  u by  v f then 

^(V0)-1^1  + V°)(v1  - v°)  = - Wp°  . 
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Thus  we  see  that  this  iteration,  which  can  be  guaranteed  to  converge  for 
the  proper  choice  of  u!  since  r is  now  a gradient  mapping,  is  quite  close 
to  the  previous  one.  The  extra  factor  in  J is  the  matrix  Hi  + (V°)-^V^)  , 
which  if  the  iterates  are  close  enough  will  be  near  the  identity. 

If  we  assume  that  the  iterates  in  (3.2)  satisfy 

V1 (V°)-1  ^ 1 + T T > 0 , 

then  the  conditions  for  using  the  approximate  relaxation  methods  of  the 
Appendix  to  this  report  are  satisfied.  In  particular, the  condition  (3.3) 
of  that  section  is sati sf i ed  wi th  3 = 2/(2  + T) . We  may  therefore  consider 
(3.2)  as  an  approximate  relaxation  method  to  (3.3).  Since  (3.2)  contains 
no  square  roots,it  may  seem  less  expensive  to  use. 

\ 
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APPENDIX  A 


AN  APPROX I MATE  BLOCK  RELAXATION 

* 

Samuel  Schechtor 


1 . I ntroduct ion 

In  the  iterative  method  of  nonlinear  relaxation  it  is  often 
desirable  to  avoid  derivative  computation  or  to  use  other  approximations 
for  the  evaluation  of  the  iterates.  In  L3]  it  was  shown  that  such 
approximations  are  valid  if  coordinate  directions  are  used  at  each  step. 
In  this  note  we  indicate  an  extension  of  this  method  to  block  relaxation 
which  would  include  a modified  Newton  method  for  a restricted  class  of 
problems . 

An  example  is  provided  to  show  that  positivity  of  the  Hessian  of 
a smooth  convex  function  does  not  imply  the  same  property  for  the  finite 
difference  approximation  in  the  large,  for  more  than  one  dimension. 


2.  Assumptions  and  Estimates 

We  will  use  the  notation  of  [2]  and  f3]  and  propose  to  minimize 

. 3 n 

real  valued  smooth  function  g(u)  ( C (R  ) by  iteration.  We  assume 

on  # o o 

u ( H is  a given  guess  to  solve  r(u)  = g (u)  = o.  Let  r = r(u  ), 

A = g"(u)  be  the  Hessian  matrix,  and  m a multiindex  taken  from  the 

set  (l,...,n)  = Z.  Let  A denote  the  principal  submatrix  of  A(u°) 

m 

defined  by  m.  We  henceforth  assume  that  A is  positive  definite. 


m 

Given  a nonsingular  matrix  K of  the  same  order  as  A , we  define 

m 

an  approximate  block  relaxation  step  by 


1 o o 
u = u + uxl  , 

(2.1)  d°  £ d = - K'V, 

m m 

o 

d ,=  o 

m 


a 
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Where  m*  is  the  complementary  multiindex  to  m and  u),  0 < U)  < 2,  is 

some  relaxation  parameter.  Since  we  wish  g(uS  to  decrease,  we 

o  o o 

require  that  (r  , d ) < 0 so  that,  for  r ^ 0, 

m 


(r  , K >•  >0. 

m . s 

We  do  not  require  K to  be  symmetric. 

To  estimate  the  change  in  g or 

1 o 

Ag  = g(u  ) - g(u  ) 

we  use  Taylor's  expansion.  Using  the  notation  B^(u)  = 3^A(u)  and 

B (u)  for  its  corresponding  principal  minor  we  get 
.1m 

/O  1 o 11  o olo 

Ag  = (g  (u  ),  u - u ) + -(u  - u , A (u  ) ( u - u )) 


1 o 


1 i u 1 O 1 O 

. Z y (u  - U , B (z)(u  -U  )(U  -U  ) 

+ 6 4 J j 


o 12  13 

= U)(r  ,d)  + - 0)  (d,  A d)  + - tu  £ (d,  B (z)d)d 
m2  m 6 jm  J 

Jem 


o 1 

for  a suitable  z between  u and  u . 


I f wo  let 


a = (d,  A d)/(d,Kd)  > 0 
m 


then 


1 12 

(2.2)  -Ag  = - ou(d,Kd)[2  - ou) U)  £ <d,B  d)d  /(d.Kd)^. 

2 3 Jem  Jm  J 

We  assume,  as  in  r3l,  that  for  W on  the  line  segment  I Joining  u°  to 

o o 
u + 2d  , 

( z |b  (w)  I2)  s ^ 

Jem  J"1 


where  the  spectral  norm  is  used  in  the  sum. 
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^ = (d,  Kd )/ (d  pd) 
K 


*K  ’ “|d|/\ 


a = max  (cr,  1) 
o 


then  i f we  set 


y ” n U 

° a + «Aor"  + - *b  ) 

o o 3 K 


it  follows  that  y s 1.  1 f we  then  choose  <u  in 

o 


0 < (D  < 2y 


wc  get  , as  in  r3"] , that 


-Ag  i - u)  (d,  Kd)  (2  - a>)  > 0 
2 ► 


where  U)  = uu/y  . This  guarantees  that  g will  decrease  with  each  iterate 
o 

whose  active  residual  is  not  zero.  To  get  the  next  iterate,  a new  m,  K and 
u)  are  to  be  chosen.  They  may  of  course  be  the  same  as  the  previous  choice 
in  certain  instances. 

3 . Convergence 

In  order  to  obtain  convergence  of  the  iterates  we  need  further 
restrictions  on  the  choice  of  the  matrices  K.  A method  for  obtaining  this 
is  to  get  an  estimate  of  the  form 


(3.1)  -Ag  2 C | r | 
m 

where  C is  a constant  Independent  of  the  Iterates.  This  requires  uniform 

upper  and  lower  bounds  on  the  K matrices  and  a lower  bound  on  y . 

o 
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This  may  be1  achieved  by  assuming  that  K Is  n reason  able  approximation 

to  A and  that  the  spectral  norm  of  K Is  uniformly  bounded.  To  this  end 
m 

we  assume  that  there  exist  positive  constants,  8 and  C independent  of 

o 

the  iterates  such  that 

(3.2)  l K | < C 

o 

(3.3)  8(d,  Ad)*  (d,  Kd). 

m 

These  inequalities  readily  yield 

2 

(3.4)  (d , Kd)  * C | r° | . 

1 m 

The  lower  bound  of  y is  obtained  from  a lower  bound  on  4>  • Since  A is 

o Km 

positive  definite  it  follows  as  in  [2]  that 

A > X >0 

m 

<J>  2 > 0 

K 

from  a priori  estimates  on  the  level  set  containing  the  iterates.  This 
is  attained  by  the  assumption  that  the  original  level  set  is  bounded. 

Under  the  further  assumptions  of  theorem  4.1  of  [3]  we  obtain  global 
convergence  of  this  approximate  block  relaxation  process.  We  may  state: 
^Theorem  3.1.  Assume  that  the  level  set 

D = [u|g(u)  < g(u°)] 

is  bounded  and  that  the  sequence  of  multi  lndeces  fm  1 is  a cyclic  ordering 

P 

covering  Z infinitely  often.  Assume  that  the  sequence  of  matrices  {a  j 

nip 

are  positive  definite  and  that  the  matrices  (K  j are  given  to  satisfy  (3.2) 

P 

and  (3.3)  then  a sequence  • ui  } may  be  chosen  so  that  for  the  process  (2.1), 
P 

r(u  )-  0.  If  the  solution  is  unique,  the  Iterates  {uP]  converge  to  the 
solution. 

We  note  that  convexity  of  g is  not  required,  only  convexity  in  the 

subspaces  defined  by  [m  ]. 

P 

The  proof,  once  the  estlaate  (3.4)  is  known  follows  along  the  lines 


of  Theorem  8.1  of  T2 "1 


4. 


Remarks 


1)  It  is  clear  from  (2.2)  that  when  the  matrices  B , Jfm  are  all 

Jm 

positive  definite  and  if  the  d < 0,  we  obtain  a decrease  in  g for 


0 < o U)  <■  2 


and  v does  not  enter.  For  the  usual  choice  of  nr  ■=  l,  the  full  range 

o 

of  u)  is  then  available.  This  was  noted  previously  T3l , in  the  scalar 
case.  An  example  of  this  is  given  by  a class  of  semilinear  elliptic 
equations  of  the  form  A<t>  = f ('f')  where  f " (<T>)  > 0.  For  the  usual  finite 
difference  approximation  to  such  problems  the  B take  the  form 

.I1" 

..  T n 

B = f (*>  ) (e  e ) , where  e is  the  jth  unit  vector  in  R,  and  are 

jm  j j j m j 

nonnegative  definite. 

o 

If  r ^ 0 and  K is  monotone  we  will  get  d < 0.  The  choice  of  K = A 
m m 

in  t tie  cited  example  will  yield  this  property. 


ii)  A common  choice  of  K,  in  general,  might  be  a finite  difference 
approximation  to  A 

m 

A .r 

a (u)  = -i-i  + 0 (h  ) 
fj  h J 


A r = r (u  + he  ) - r (u). 

Ji  1 j J i 


Thus 


K = A + BH  = (A  r /h  ) 

■n  J i J 


where  H is  a diagonal  matrix  with  entries  h , jfm  and  B contains  entries 
in  g'''.  If  H is  sufficiently  small,  we  may  estimate 


|(d,  BHd ) | S f(d,d) 


where  f depends  only  on  g.  We  then  get 
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(d,Kd)  S (d,  A d)  (1  - */<&  ) 
m A 

>0  (d,  Ad) 

m 

with  0=1-  f /X , which  is  positive  for  ( < \. 

ill)  If  we  wish  to  choose  K = I,  the  identity,  then  we  may 

choose  6 = 1/A  where  A is  the  largest  of  the  eigenvalues  of  A in 

m 

a suitably  large  but  bounded  domain. 

iv)  Since  the  method  described  includes  the  case  of  m being 
the  full  index  set  (1,  2,...,n),  we  may  consider  solving  an  elliptic 
system 

r(u)=Lu-f=0 
by  a sequence  of  problems 

K(u*  - u°)  = - u)r(u°) 

and  the  above  results  will  apply.  For  example  the  choice  of  K = - A , 

h 

the  discrete  Laplacian,  has  been  proposed  by  Concus  and  Golub  [l]  for 
solving  a linear  elliptic  problem.  This  has  the  virtue  that  K may  be 
a much  simpler  operator  to  use  than  L.  Our  results  indicate  that  such 
methods  are  feasible  in  the  nonlinear  case  as  well. 

5.  A Counterexample 

For  a smooth  convex  function  of  one  variable,  it  is  well  known 
that  the  difference  quotient  of  its  derivative  is  nonnegative: 

(g'(u+h)  - g’(u))/h  S 0 

for  all  u and  u + h in  the  domain  of  definition.  Furthermore  if  g"(u)  & \ > 0 
for  all  u,  then  its  difference  quotients  have  the  same  lower  bound  for  all 
u and  h. 

The  question  now  arises  about  the  corresponding  statement  in  higher 

ft 

dimensions.  That  is,  if  we  define  the  matrix 
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T 


D . <4lJ) 

by  A = A r /h  can  we  conclude  that  (w,  I>w)  > 0 for  all  u and  h , w 

ij  j i .J  J 


If  A(u)  ^ 0 for  all  u? 

Furthermore  if  A(u)  > \ > 0 globally  can  we  conclude  that  there 

is  a X > 0 such  that 
o 

(w,  Dw)^  X (w,w) 

o 

for  all  w,  h^  ? 

It  follows  from  the  previous  section  that  for  sufficiently 
smooth  g the  statements  are  true  for  all  u but  for  sufficiently 
small  h . We  will  show,  by  a counterexample,  that  the  statements 

j 

are  not  valid  for  nil  h.. 

J 2 2 $ 

The  example  we  will  first  use  is  g = (x  + y ) in  two 

dimensions.  Thus  g is  convex  and  r = x/g,  = y/g.  We  choose 

the  point  (x  , v ) = (1,  1)  and  set  h = h,  h = 0 so  that  A,_  = a , 
o o 12  12  12 


A 

o o 

a 

22 

22 

/A 

a 1 

/ 11 

12\ 

D =| 

\A 

a ) 

\ 21 

22/ 

K = 

g(l.l)  = -h  , 

e.  = i 

o 

1 

A 

<!>  - , 

A 

11 

g,  g 

21 

°12 

1 o 

-3 

- *yg  , 

a = 

22 

If  it  were  true  that 

of  D 

would  be  nonnegative 

that 

det  (D)  *■*  0 and 

get 

a = 

-3 

g = - a 

22 

o 12 

2 -3 

{ R . 


i 
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v 1 , 1 2+h  , 

dct(D)  = — (A  + A > = — ( .72) 

J 11  21  .36, 

K hg  1 

° 2 ° 

2 (1  4 x ) x 

2+h  1 1 

But  ( ) = — rr  1 4 2 < 2 

g 2 2 

*1  1 4 X 1 4 X 


for  h ^ 0,  so  that  the  det  (D)  < 0 and  we  get  a contradiction. 

00 

To  obtain  an  example  for  a smooth,  in  fact  C , convex  function 
we  choose  a perterbation  of  g: 

i 

2 2 J 

g = (x  4 y 4 f)  for  e > 0. 

( 

Let  D be  the  corresponding  matrix  at  (1,1).  It  follows  from  continuity 

e 

considerations  that  the  det(D  ) is  negative  for  small  enough  f.  In  fact 

e 

a direct  check  for  ( S 1/2  and  h > 20f  shows  that  det(D  ) < 0.  Since  the 

e 

Hessian  of  g is  uniformly  bounded  from  below,  both  coniectures  are  false. 

€ 

It  follows  from  the  remark  (ii)  that  if  we 
is  positive  definite  and  of  order  |hb|,  we  will 


choose  K = D 4 E,  where  E 
satisfy  the  positivity. 
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PUBLICATIONS  SUPPORTED  BY  THIS  CONTRACT 

Schechter,  S.,  "On  the*  Choice  of  Relaxation  Parameters  for  Nonlinear 
Problems"  in:  Numerical  Solution  of  Systems  of  Nonlinear  Algebraic 
Equations,  eds.  G.  D.  Byrne  anti  C.  A.  Hall,  Academic  Press,  New  York, 
pp.  3-49-372,  1973. 

Schechter  S.,  "An  Approximate  Block  Relaxation"  to  appear  in  SIAM 
J.  Numer.  Anal. 

Ablow,  C.  M.  and  S.  Schechter,  "Compylotropic  Coordinates",  submitted 
for  publication,  Abstract  in  Notices  of  the  American  Mathematical 
Society,  Vol . 24,  No.  3,  No.  746-C3 , 1977. 
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